library(R.matlab)
library(ggplot2)
library(lme4)
# library(lmerTest)
library(nlme)
library(multcomp)
library(car)
library(ez)

# Set paths based on current machine
compname = Sys.info()
compname = compname[[4]]
if (compname == "DESKTOP-BHR0AU7") {
  loadPath = 'E:/ANL/Experiments/RESULTS/VAST/DATAforR/'
} else if (compname == "MSI") {
  loadPath = 'E:/ANL/RESULTS/VAST/DATAforR/'
} else {
  stop('set directories for this machine')
}

# Load and format the ERP data
erp_mat = readMat(paste(loadPath, 'ERP_peaks.mat', sep=""))
ERP_peaks = erp_mat$ERP.peaks
condnames_INT = unlist(erp_mat$cnames.int)
Nss = nrow(ERP_peaks)
WM_mod = unlist(lapply(condnames_INT, substr, 1, 1))
WM_dom = unlist(lapply(condnames_INT, substr, 2, 2))
int_task = unlist(lapply(condnames_INT, substr, 4, 5))

ERP_peak_df = data.frame('subj'=rep(1:Nss, length(condnames_INT)), 'mod'=rep(WM_mod, each=Nss),
                         'dom'=rep(WM_dom, each=Nss), 'int'=rep(int_task, each=Nss),
                         'peak'=c(ERP_peaks))
ERP_peak_df$subj = factor(ERP_peak_df$subj)
ERP_peak_df$mod = factor(ERP_peak_df$mod)
ERP_peak_df$dom = factor(ERP_peak_df$dom)
ERP_peak_df$dom = factor(ERP_peak_df$dom, levels=levels(ERP_peak_df$dom)[c(2,1)])
ERP_peak_df$int = factor(ERP_peak_df$int)
ERP_peak_df$int = factor(ERP_peak_df$int, levels=levels(ERP_peak_df$int)[c(2,1)])


# Construct the models
mdl = lmer(peak ~ mod * dom * int + (1|subj), data = ERP_peak_df,
          control=lmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e6), calc.derivs=FALSE))
mdl_no3way = lmer(peak ~ mod + dom + int + mod:dom + dom:int + (1|subj), data = ERP_peak_df,
           control=lmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e6), calc.derivs=FALSE))
# Assess significance of the fixed effects with an anova
mdl_stat = Anova(mdl, type='III')

## --- Helper function for generating fit check plots --- ##
check_model_fit = function(model)
{
  windows()
  print(plot(model))
  windows()
  qqnorm(resid(model))
  qqline(resid(model))
}


# Basic 3-way ANOVA version
ERP_anova = ezANOVA(ERP_peak_df, dv=peak, wid=subj, within=.(mod, dom, int))

# Out of curiosity, what does it look like with just the visual WM task data?
ERP_df_visOnly = ERP_peak_df[ERP_peak_df$mod=='v',]
ERP_df_visOnly = ERP_df_visOnly[, names(ERP_df_visOnly) != 'mod']
ERP_anova_visOnly = ezANOVA(ERP_df_visOnly, dv=peak, wid=subj, within=.(dom, int))
# Tukey's posthoc testing
ERP_df_visOnly$cInter = interaction(ERP_df_visOnly$dom, ERP_df_visOnly$int)
visOnly_lme = lme(peak ~ cInter, random = ~1|subj, correlation=corCompSymm(form=~1|subj), data=ERP_df_visOnly)
visOnly_posthoc = summary(glht(visOnly_lme, linfct=mcp(cInter = "Tukey")), test = adjusted(type = "holm"))
print('P2 Post-Hoc Tests (Visual WM conditions only)')
print(visOnly_posthoc)

# And only auditory data
ERP_df_audOnly = ERP_peak_df[ERP_peak_df$mod=='a',]
ERP_df_audOnly = ERP_df_audOnly[, names(ERP_df_audOnly) != 'mod']
ERP_anova_audOnly = ezANOVA(ERP_df_audOnly, dv=peak, wid=subj, within=.(dom, int))
# Tukey's posthoc testing
ERP_df_audOnly$cInter = interaction(ERP_df_audOnly$dom, ERP_df_audOnly$int)
audOnly_lme = lme(peak ~ cInter, random = ~1|subj, correlation=corCompSymm(form=~1|subj), data=ERP_df_audOnly)
audOnly_posthoc = summary(glht(audOnly_lme, linfct=mcp(cInter = "Tukey")), test = adjusted(type = "holm"))
print('P2 Post-Hoc Tests (Auditory WM conditions only)')
print(audOnly_posthoc)



                            